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Abstract 

Collagen fibres play an important role in the mechanical behaviour of many soft tis- 
sues. Modelling of such tissues now often incorporates a collagen fibre distribution. How- 
ever, the availability of accurate structural data has so far lagged behind the progress 
of anisotropic constitutive modelling. Here, an automated process is developed to iden- 
tify the orientation of collagen fibres using inexpensive and relatively simple techniques. 
The method uses established histological techniques and an algorithm implemented in the 
MATLAB image processing toolbox. It takes an average of 15 seconds to evaluate one 
image, compared to several hours if assessed visually. The technique was applied to his- 
tological sections of human skin with different Langer line orientations and a definite cor- 
relation between the orientation of Langer lines and the preferred orientation of collagen 
fibres in the dermis (P<0.001, R^= 0.95) was observed. The structural parameters of the 
Gasser-Ogden-Holzapfel (GOH) model were all successfully evaluated. The mean disper- 
sion factor for the dermis was K = 0.1404 ± 0.0028. The constitutive parameters /i, ki and 
^2 were evaluated through physically -based, least squares curve-fitting of experimental test 
data. The values found for /i, ki and ^2 were 0.2014 MPa, 243.6 and 0.1327, respectively. 
Finally, the above model was implemented in ABAQUS/ Standard and a finite element 
(FE) computation was performed of uniaxial extension tests on human skin. It is expected 
that the results of this study will assist those wishing to model skin, and that the algorithm 
described will be of benefit to those who wish to evaluate the collagen dispersion of other 
soft tissues. 
Keywords: Fibre Orientation, Anisotropic, Skin, Collagen Fibres 



1 Introduction 

Collagen fibres govern many of the mechanical properties of soft tissues, in particular their 
anisotropic behaviour. Due to the complex nature of fibre arrangements, such tissues are often 
represented as either isotropic or transversely isotropic [[6l[2l[3l. The availability of quantitative 
structural data on the orientation and concentration of collagen fibres is crucial in order to 
describe the behaviour of these soft tissues accurately. 

A number of studies have used statistical distributions to describe the fibre arrangement in 
soft tissues. Lanir [20J was the first to attempt to account for fibre dispersion. Lanir's method 
expresses the mechanical response in terms of angular integrals. This technique accounts for 
the contribution of infinitesimal fractions of collagen fibres in a particular orientation. This ap- 
proach leads to accurate results but it is not a practical option for efficient numerical implemen- 
tation JH. The second approach uses generalised structure tensors (GST), which are assumed 
to represent the three-dimensional distribution of the fibres. The strain energy function is then 
calculated by using the average stretch rather than by using a stretch for each individual fibre. 
GST is a simpler method than Lanir's, and is easily implemented in FE algorithms [ [T0) |. 

In this paper we implement the GST method proposed by Gasser et al [lOJ, described in 
detail in Section [23| In that model, structural parameters such as the fibre dispersion parameter 



and mean orientation of fibres are required. In order to quantify the orientations of collagen fi- 
bres in the human dermis an automated process is developed for detecting collagen orientations 
in histology slides. 

Histology is the microscopic study of cells and tissue. It is an important diagnostic tool in 
medicine which is utilised here for the purpose of analysing collagen orientation. Traditionally, 
histology slides are examined by expert observers who assess individual slides visually. This 
method is both time consuming and subjective. An automated process would enable large vol- 
umes of images to be analysed quickly and improve the objectivity of the task. Much research 
has been conducted recently on imaging techniques of biological soft tissue for the extraction 
of structural data; however these techniques require either expensive equipment or manual post 
processing [[8l|30l|3ll|32]|. There has been little research to date carried out on the automated 
analysis of histology slides. Four notable exceptions are the work of Van Zuijlen et al ll29l . 



Noorlander et al Jll], Elbischger et al (5] and Jor et al ifTTl . 

Elbischger et al [5 J successfully automated the analysis of collagen fibre orientation in the 
human adventitia, (the outermost membrane of the artery). Their technique was based on the 
automated analysis of Transmitted Light Microscopy images which were stained with Elastica 
van Gieson. The algorithm uses a ridge and valley analysis to detect the orientation of collagen 
fibres and segments regions of homogeneous fibre orientations. This is a somewhat complicated 
technique which requires substantial coding to implement. There are also assumptions made 
which cannot be applied to the dermis, i.e. that the image contains at least 50% collagen and 
that the fibres have a common orientation. In contrast, Noorlander et al f24l used a relatively 
simple technique upon which our own algorithm is based. In their study, Picrosirius red staining 
was used together with epipolarised light to image the fibres. Individual collagen fibres were 
detected and ellipses fit to the longest ten fibres. While their technique is quantitative it is not 
an automated technique and requires the user to identify fibre orientations visually. The results 
refer to the 'collagen alignment index' which is the mean length of the best fit ellipse and has 
no connection with the angular orientation of fibres, which itself is not measured. Van Zuijlen 
et al ll29ll used Fourier analysis to measure the level of anisotropy of collagen in the skin. 
Those authors concluded that analysis of 'orientation index' by Fourier analysis is superior 
to conventional techniques. However, the 'orientation index' is a measure of the anisotropy 
of the matrix but, just like the 'collagen alignment index' of Noorlander et al ll24l . it fails to 
provide information on the mean orientation of the fibres. In a recent publication Jor et al 
ifTTIl modelled the orientation of collagen fibres in porcine skin. Their technique involved the 
staining of porcine skin sections with Picrosirius red and image acquisition by confocal laser 
scanning microscopy. Their study provided quantitative results on the orientation of collagen 
fibres in porcine skin, however the plane of interest was normal to the epidermis and it is 
generally believed that the preferred orientation of collagen fibres is parallel to the epidermis 

The chief advantage of the technique proposed here over other automated techniques is that 
it is an inexpensive and relatively simple technique. It can be easily implemented in MATLAB 
using the Image Processing Toolbox so that the algorithm can be easily amended as required 



for the user's specific application. It is a fully automated process (aside from the image ac- 
quisition phase) and is capable of analysing a greater number of images, and faster, than a 
manual analysis. It also eliminates the subjectivity which is present using traditional meth- 
ods. The technique provides quantitative data on both the collagen orientation and the level of 
anisotropy in the dermis, but it could be easily amended to identify any tissues/objects within 
other soft tissues. 

Due to the anisotropic properties of collagen based soft tissues ll23]| . modelling now often 
incorporates a collagen fibre distribution. The accuracy of these structural models rely heavily 
on knowledge of collagen fibre orientations. While this data has recently been published for 
porcine skin by Jor et al [fTTl . to the best of the authors' knowledge, quantitative data on the 
orientation of collagen fibres in the human dermis has not been previously published. With the 
quantitative data obtained here, structural parameters such as 7 and K from the Gasser-Ogden- 
Holzapfel model ifTOl can be evaluated independently. This structural data, coupled with the 
mechanical testing of the same skin samples ll23]| provides us with sufficient data to model 
human skin using the GOH model. 

2 Materials & Methods 
2.1 Tensile tests and histology 

In vitro tensile tests of human skin were performed in Ifsttar (Institut frangais des sciences 
et technologies des transports, de I'amenagement et des reseaux), France. French law allows 
human corpses that have been donated to science to be used for research purposes. The ethics 
committee within Ifsttar approved the use of human biological material. 

The tensile tests were performed using a Universal Tensile Test machine at a strain rate 
of 0.012s~^ The tensile load was measured with a IkN piezoelectric load cell and the strain 
was measured via a displacement actuator. Skin biopsies were excised adjacent to the tensile 
test samples and stained with Van Gieson to dye collagen fibres pink/red. Further details of 
the tensile test experimental procedure and histology protocol are outlined in Ni Annaidh et al 
lIZl. 



2.2 Automated detection of fibre orientations 



Images were taken of slides parallel to the epidermis (see Section [2!4| ) using an Aperop ScanScope 



XT scanner and ScanScope software. This slide scanner has the ability to scan up to 120 slides 
automatically. The images to be analysed were selected from the reticular layer of the der- 
mis which forms the main structural body of the skin. The optical magnification used was 5x, 
which is quite low; however, this was the optimal magnification to obtain a large enough field 
of view to be representative while also maintaining enough detail to recognise the boundaries 
of fibre bundles. Multiple images were taken of each sample in order to capture the entire area. 
The field of view for each image taken at this magnification was 2.24 mm^. The orientation of 
collagen fibres were then calculated in a fully-automated customised MATLAB routine using 
the Image Processing Toolbox. The algorithm is described below and is illustrated in Fig.[T] 

• A global threshold level was computed using the graythresh function. This function 
chooses a global image threshold using Otsu's method which aims to minimise the intr- 
aclass variance of black and white pixels ll26ll . This automatically distinguishes collagen 
fibres from other areas which are not of interest in the analysis, such as cells or histo- 
logical ground substance. The image was then binarised based on this threshold level 
producing an image containing black pixels for collagen and white pixels for all other 
areas; 

• Morphological operations were performed on the binary image using the bwmorph func- 
tion. An erosion step was performed to detach cross-linking fibres from each other. This 
step removed pixels from the boundaries using a structural element of size 3 pixels x 
3 pixels. One iteration only of this step was performed which was the optimal number 
of iterations to detach cross-linking fibres while also leaving smaller fibres intact. The 
second morphological operation was performed using the imfill function. This function 
fills in the 'holes' within the binary image, where a hole is defined as a set of isolated 
pixels which cannot be reached by filling in the background pixels from the edge. This 
step resulted in a 'cleaner' image to analyse. 

• Individual fibre bundles are identified using bwlabel. This function identifies connected 



components (8-connected) in a 2D binary image and labels each component individually. 
It uses the general procedure outlined by Haralick and Shapiro [fT2||: 

• The regionprops function outputs a set of properties for each labelled component. The 
function also fits an ellipse to each component by matching the second order moments of 
that component to an equivalent ellipse following the procedure by Haralick and Shapiro 

El; 

• Only ellipses which were elongated and of a certain area were selected Q 

• The orientation of the major axis of each ellipse was calculated and taken as the ap- 
proximate orientation of each component. The advantage of fitting an ellipse about each 
component is that the orientation is measured in a systematic and repeatable manner 
which can account for the non-uniformities of the shape of each component; 

• The orientations of each component was then plotted on a histogram (see Fig. [2]). Two 
distinct peaks were evident from this figure. It is assumed that these two peaks corre- 
spond to the preferred orientation of two crossing families of fibres as shown in Fig.|3j A 
von Mises probability density function was fit to the data and the mean orientation and 



dispersion factor were calculated as described in Section 2.4 



Validation 

As explained by Van Zuijlen et al ll29l , there are difficulties involved with validating automated 
techniques because the true result of each histology slide is unknown. For the validation of their 
algorithm, computer generated images representing collagen fibres with known orientations 
were used. However, in order to create an accurate computer generated image, a model which 
represents the structure of the collagen matrix is required, but no such model exists JSl. To 
validate our technique, a selection of slides in the plane perpendicular to the epidermis (see 



^Out-of-plane fibres manifest themselves as circular areas. This condition was introduced to exclude out-of- 
plane fibres from the analysis. A sensitivity analysis was carried out to determine the sensitivity of the algorithm 
to varying the area and eccentricity criteria. It was found that a ±10% change in both the area and eccentricity 
criteria led to a ±1% change in the mean orientation. This indicates that the selected criteria do not significantly 
affect the result of mean orientation. For this study it was specified that in-plane fibres must have an eccentricity 
larger than 0.7 and an area greater than 1000 pixels, which for our images, captured at 5x, corresponds to 5^m^. 



Fig. |4]) were manually segmented and their mean orientation compared to those calculated 
through the automated process. Although it is the plane parallel to the epidermis that is used 
for the collection of structural data, the perpendicular plane was chosen for validation purposes 
because in general it displays a more orientated pattern and is easier to segment manually. 
Manual segmentation was performed by marking bundles of elongated fibres with a line as 
shown in Fig. [5} The orientation of each line was then measured manually and the mean 
calculated. 

The dermis possesses a complex interwoven collagen structure, which is difficult to assess. 
It is quite possible that the collagen pattern may change over a small volume. For this reason 
six images were taken at each depth, whereby each image spanned 3J5mm^ of the SOmrn^ 
sample. The analysis was performed at six levels of increasing depth and the mean orientation 
was taken as the average over these six levels. 

2.3 Gasser-Ogden-Holzapfel Model 

The Gasser-Ogden-Holzapfel (GOH) model applies to incompressible solids with two preferred 
directions aligned along the unit vectors a\ and a2 (say) in the reference configuration (see 
Fig.|3]). Its strain energy density ^ is of the form 

^ = ^(C,Hi,H2) (1) 

where C is the right Cauchy-Green strain tensor, and the structure tensors Hi , H2 depend on 
ai and a2 and on the dispersion factors K*i, K2 (to be detailed later), respectively, as follows 

H/ = K-/I+(l-3K-/)ai®a2, (/= 1,2). (2) 

Specifically, the GOH model assumes that ^ depends on /i = tr(C), tr(HiC) and tr(H2C) 
only, as follows 

^=^(/i-3) + M I ^|e^-KH,c)-ip_i\ (3) 



where ji, kn, fc/2 are positive material constants, and from Eq.|2| 

tr(H,-C) = K-,-/,- + (l-3/cO/4/, (4) 

with l4i = a/Ca/, two anisotropic invariants. Note that the constitutive parameter ii has the 
dimensions of stress: it would be the shear modulus of the solid if there were no fibres (kn = 0); 
whilst the parameters kn and ki2 are dimensionless stiffness parameters: the kn are related to 
the relative stiffness of the fibres in the small strain regime, and the ki2 are related to the large 
strain stiffening behaviour of the fibers. 

Now we focus on homogeneous uniaxial tensile tests. These can be achieved for anisotropic 
tissues when two families of fibres are mechanically equivalent kn = /:2i = h (say) and ku = 
kii = ^2 (say), with the same dispersion factor K*i = K*2 = K* (say), and when the tension occurs 
along the bisector of ai and a2 (see Fig. [3]). Let us call 7 the angle between ai and the tensile 
direction, so that now 

ai = cosyi + sinyj, a2 = cosyi — sinyj. (5) 

Here, i is the unit vector in the direction of tension, and j is the unit vector in the lateral 
direction, in the plane of the sample. The stretch ratios along those unit vectors are Ai and A2, 
respectively. Then /41 = I42 = ^f cos^ 7+^1 ^^^^ Y = h (^^j) ^^^ ^ reduces to 

.p = ^(/j _3) +^^ |e^2['c/,+(i-3^)/4-i]2 _ jl (6) 

2 k2 ^ J 

giving the following expression for CJ, the Cauchy stress tensor 

a = -pl + 2--FF^ + -- [Fai®Fai + Fa2®Fa2] , (7) 

oil 0I4 

where /? is a Lagrange multiplier introduced by the internal constraint of incompressibility and 
F is the deformation gradient. Note that Fai®Fai +Fa2(8)Fa2 = 2(Ai cos y)H 1 + 2(^2 sin 7)^j (g) j. 



showing that a is diagonal in the {i, j, k} basis. Its components are 



Oil = -p + 2(^1 + ^4 sin^ r)^2 = 0, 

0-33 = -p + 2^1 Af%-^ = 0, (8) 



where 



.^2«^^ 



2^i=A^(l+4fcime^2« )^ 

2^4=4M^i(l-3/c)ae^2«', 

a = /c(Af + A| + Af%-^) + (1 -3/c)(Af cos^7+A|sin^7) - 1. (9) 

Now, ehminate p from the stress components to get the two equations 

cjii =A^(Af-Af%-^) +4^1)^1 ae^^a' [/c(Af - Af%-^) + (l -3/c)Af cos^y] , (10) 
= A| - Af %-2 + 4ytiae^2a^ [/c(A| - Af %-2) + (1 - 3/c)A|sin2 7] . (11) 

Equation ( [TT] ) gives the relationship between the tensile stretch and the lateral stretch and 



allows, implicitly, A2 to be expressed in terms of Ai. Substituting then into ( ^0| ) gives the CJi 



Ai stress-stretch relationship. In the isotropic limit, K* = 1/3 (see Section 2.4), and ( [TT] ) yields 

— 1/2 

the well-known relationship A2 = A. ^ for uniaxial tension in incompressible solids. 



The two equations ([TO|)-([TT]) form the basis of a numerical determination of the constitutive 



parameters /i, ki and k2, assuming that the structural parameters K and 7 are known. It should 
be noted here that the inclusion of fi in the anisotropic term of equation ([6]) is not standard for 
the GOH model and has been added here for ease of calculation. We now quantify further those 
latter parameters, K and 7 . 
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2.4 Fibre Dispersion 

The GOH model assumes that the mean orientation of collagen fibres has no out-of-plane com- 
ponent. Our histological examination of the skin indicates that the majority of collagen fibres 
in the dermis run parallel to the epidermis. Slides parallel to the epidermis had, on average, 
three times less cross-sectioned fibres than slides perpendicular to the epidermis. This is in 
agreement with Holzapfel et al lfT3]| who state that the preferred orientation of the 3D collagen 
fiber network lies parallel to the surface, but to prevent out-of-plane shearing, some fiber ori- 
entations have components which are out-of-plane. Despite the assumption that the fibres have 
no out-of-plane component in the GOH model, the three-dimensional nature of the adopted 
distribution implies that although the preferred orientation of the fibers are in the plane parallel 
to the epidermis, some fibers orientations have an out-of-plane component ifTOll . 

Here we assumed that each of the two families of collagen fibres is distributed according 
to a ;r-periodic Von Mises distribution, which is commonly assumed for directional data. The 
standard ;r-periodic Von Mises Distribution is normalized and the resulting density function, 
p(0), reads as follows, 

P(e)^47A --p"";°;'g'^" . (12, 

V 271 erfi(v2Z?) 
where b is the concentration parameter associated with the Von Mises distribution and is the 
mean orientation of fibres (for a graph of the variation of the dispersion parameter K with the 
concentration parameter b, see Gasser et al lITOl ). 

The parameters b and were evaluated using the rule function in MATLAB. Analogous to 
least squares curve-fitting, the maximum likelihood estimates (MLE), is the preferred technique 
for parameter estimation in statistics. Then K is calculated by numerical integration of the 
integral given by Gasser et al ifTOll . 

1 r^ 
K-=- / p(0)sin^0J0. (13) 

4 Jo 

The structural parameter K describes the material's degree of anisotropy. It must be in the range 
^ K* ^ 1/3: the lower limit, K* = 0, relates to the ideal alignment of collagen fibres and the 
upper limit, K* = 1/3, relates to the isotropic distribution of collagen fibres. Fig. |6]is a 3D 
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graphical representation of the orientation of collagen fibres for different values of K. 

The fibres were assumed to form an interweaving lattice structure as first postulated by 
Ridge and Wright fTJ] and shown in Fig. |3| These authors suggested that the mean angle of 
the two families of fibres indicates the direction of the Langer lines. More recent in vitro ifTTl 
and in vivo ll28]| studies have also supported this hypothesis. The lattice structure proposed by 
Ridge and Wright ETl is an idealised one, and the adoption here of the dispersion factor creates 
a more realistic scenario. Finally, recall that the two families (directions) of fibres are assumed 
to have a common dispersion factor. 

2.5 Finite element representation 

An FE computation was used to simulate uniaxial tensile tests of human skin which were 
described in a previous publication Il23ll . The simulation was carried out for three samples: 
parallel, perpendicular, and at 45° to the Langer lines. The test samples were of the dimen- 
sions shown in Fig. [7| The length of the undamped specimen was 68 mm and the thickness 
was 2.25 mm. 1512 reduced integration hybrid hexahedral (C3D8RH) elements were used 
for the mesh. The numerical analyses were performed using the static analysis procedure in 
ABAQUS/Standard. Displacement was applied through a smooth amplitude boundary condi- 
tion, and the top and bottom of the sample were encastred to represent the clamping of samples. 
The material model used was the anisotropic GOH model, which is an internal material model 
inABAQUS. 

3 Results 

3.1 Structural parameters 

As expected, it was found from the histology that the collagen fibres were locally orientated. 
This meant that each sample had a different mean orientation and fibre dispersion. Table [T] 
tabulates the results of 12 different human skin samples, with different orientations, which 
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were procured from the backs of two different subjectq^ See Fig. [s] for details of specimen 
location and orientation. The preferred orientation, 0, refers to the bisector of the two families 
of fibres. The 95% Confidence Interval of the mean of these images ranged from 1.16° — 
2.77°, thereby indicating that the preferred orientation does not change significantly over this 
small area. The results of our validation revealed that the automated process differed from the 
manual segmentation by an average of 5° ±4°, giving us further confidence in the validity of 
our technique. 

A Pearson correlation test was carried out to test for a correlation between the measured 
preferred orientation obtained through histology and the perceived orientation of Langer lines 
(the natural lines of tension in the skin). The orientation of Langer lines was assessed using 
generic maps, described further in Ni Annaidh et al ll23]| .The correlation was deemed to be 
significant (P<0.001) with an R^ value of 0.95. This shows that the Langer lines have an 
anatomical basis, a point which had previously been suggested but until now had not been 
quantitatively assessed. 

3.2 Constitutive parameters 

It was assumed that the orientation of collagen fibres is symmetric about the axis of applied 
stress. In reality, this is not always the case; however some assumptions must be made in 
order to ensure that the constitutive relations remain practical for numerical implementation. 
In particular, as explained in Section [23} this assumption leads to a homogeneous deformation 



of the sample, and in turn, to an explicit stress-strain solution. In this section, three illustrative 
examples from Table [T] (highlighted) have been chosen for further investigation. Two of these 
examples have been chosen because they are the samples that are closest to being symmetrical 
about the axis of applied stress. The third sample has been chosen to illustrate how its behaviour 
can be modelled using FE analysis and is described further in Section [331 

The constitutive parameters for the GOH model are obtained by using Equations ( [T0| ) and 
( [TT] ), obtained in Section 23 When linearized in the neighbourhood of small strains. A/ c^ 
^It should be noted here that these results have been obtained using the mle method described in Section 



i |2.4 

An alternative, but less reliable method exists whereby data is clustered into two intersecting fibre families using 
an agglomerative clustering algorithm such as was performed in Ni Annaidh et al [i23il . 
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1 + £/, say, we find that they read as follows, 

cjii =4^[l+2^i(l-3K-)^cos'^7]£i+2^[l+4fci(l-3K-)^sin^7cos^7]82, (14) 

0=[l+4iti(l-3/c)^cos^7sin^7]8i+2[l+2iti(l-3K:)^sin4 7]£2. (15) 

These expressions reveal that the constitutive parameters ji and ki are related to the early stages 
of the tensile tests, whilst ^2 is a stiffening parameter, related to the latter (nonlinear) stages of 



the tensile tests. By solving ( [15] ) for £2, and substitution into ([14]), we find the linear stress-strain 
relation Gn=E\£\, where Ei is the infinitesimal Young modulus in the 1 -direction, found here 
as 



3 + 8^1 (1 - 3k:)^(1 - 3cos^ 7sin^ 7) 
l+2ki{l-3K)^sin^Y 



El = ;, ,. / \^.., . 4.. — -1^^ (16) 



(which is consistent with the formula £" = 3/i in linear isotropic (k* = 1/3) incompressible 
elasticity). Hence, by plotting the values of cJn for the early part of the tests (first 1000 data 
say, corresponding to a tensile stretch of less than 2%), we can determine Ei by linear regression 
analysis, see Figj9fa). Here we have plotted the 'parallel' sample highlighted in Table [T] data 
for which was collected from tensile tests of human skin samples |[23]| . 

Once El is determined, /i can be expressed in terms of £"1 and ki using Eq. ( [T6] ). Then the 
remaining material parameters ki and k2 are found through the nonlinear least squares fitting 
with experimental test data of Equation ( [TO] ), subject to the definition of A2 in terms of Ai given 
by Equation ( [TT] ). The data fitting was performed using the Isqnonlin MATLAB routine in the 
Optimisation Toolbox where the objective function, Err{k), was given as 

Err{k) = Y^{yfP-y'^''''^^^f (17) 



/=i 



Where n is the number of experimental data points, y- ^ is the experimental value and 
^mo ^ I J j^ ^j^^ value predicted by the model using the current material parameters, k. 

Non-linear optimisation procedures are often sensitive to the initial starting point provided 
by the user [[25]| . In our case, the initial estimate for k\ was found by calculating the slope 
of the non-linear part of the stress-stretch curve. Since we have shown that k\ is related to 
the stiffening stage of the tensile test, our initial estimate therefore has a physical meaning. 

14 



Furthermore, the initial estimates of both k\ and ^2 were varied over a large range (0-le^ for k\ 
and 0-le^ for k2) and lead to the same set of optimal parameters each time, illustrating that the 
results of the optimisation procedure are not sensitive to this initial estimate. 

The results of our optimisation procedure gave us a value of 243.6 for k\ and 0.1327 for 
^2, with an R^ of 99.5%. FigJ9[b) shows the GOH model fit to the experimental data. It can 
be seen that these material parameters provide an excellent fitting to the 'parallel' sample, at 
least from a descriptive point of view. Turning now to the predictive capabilities of the GOH 
model, we examine a sample 'perpendicular' to the Langer lines. In Fig. [10] we use the material 
parameters k\ and k2 obtained through least squares fitting, ii calculated by linear regression 
and Eq. ([16]), coupled with the unique structural data for the 'perpendicular' sample in Table |2| 
We compare the model prediction for a tensile test occurring perpendicular to the Langer lines 
to the experimental data. The fit remains good with R^ = 97.96%. This shows that the model is 
capable of predicting the behaviour of skin once the structural parameters have been evaluated. 

3.3 Finite element simulations 

The conventions used by ABAQUS are related to ours through C\o = /i/2, k[ =k\ji and ^2 = ^2, 
so that the ABAQUS parameters used were Cio = 0.1007 MPa, k[ = 24.53 MPa and k'2 = 
0.1327. 

The FE simulation results of the uniaxial tensile tests for both the parallel and perpendicular 
samples were identical to the analytical solution. The results were independent of both mesh 
density and element type. Fig. [TT] shows the Cauchy stress distribution across the parallel and 
perpendicular samples at the end of the test. The large difference in magnitudes between the 
two is due to the variation in the mean orientation of fibres. Note the uniform distribution of 
stress in the middle section of the test samples, thanks to the dog-bone shape of the specimen, 
and the symmetry of fibres about the axis of applied stress. 



As discussed in Section [X2l for ease of determining the material parameters, it was assumed 
that the orientation of collagen fibres is symmetric about the axis of applied stress. This is 
for an idealised scenario only, where one knows the exact orientation of collagen fibres prior 
to testing, and can therefore apply the stress in this orientation. However, with the model 
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parameters that have now been determined, a FE simulation can provide results for samples 
where the collagen fibres are not symmetric about the axis of applied stress. Examining the non- 
symmetric example in Fig.[T2f a), we can see a non-uniform distribution of stress throughout the 
test specimen. A local magnification of stress occurs near the neck regions of the test specimen. 
This is due to the non-symmetry of collagen fibres about the axis of applied stress and makes 
this problem a much more complicated one to solve analytically. The presence of significant 
levels of shear in Fig.[T2]^b) (which is absent from both the parallel and perpendicular samples) 
indicates further the effect of this non-symmetry on the sample response. 

Because the stress distribution in this sample is non-uniform, to examine the predictive 
capabilities of the GOH model here, we must plot the experimental force-displacement data 



against the values predicted by ABAQUS for a node at the top of the test sample, see Fig. [13 
Again, we have found that the model predicts the behaviour well with an 7?^ of 94.4%, showing 
that the GOH model is capable of predicting the anisotropic response of human skin. 

4 Discussion 

For this paper, histology slides in three different planes were examined; however, after captur- 
ing the images from all three planes it was observed that three times as many cross-sectioned 
fibres were present in the plane normal to the epidermis, therefore an assumption was made 
to ignore the fibres normal to the epidermis. Hence, the further analysis of the samples was 
restricted to the plane parallel to the epidermis alone, making this a 2D analysis. Physically 
however, there are a percentage of fibres that run normal to the epidermis and this information 
has not been captured here. Furthermore, unloaded collagen fibres have a crimped nature: Their 
relative orientation may vary with respect to the epidermal plane. Here we have tried to over- 
come this limitation by taking an average measure over six levels of depth spanning 30/im,with 
the view that the general orientation of the collagen fibres are still captured. Ideally, a full 3D 
analysis of the dermal structure could consider the effects of collagen crimping. A 3D analysis 
was beyond the scope of this paper, but the current technique could be extended by creating a 
montage of overlapping images through the thickness of the dermis, as described in Jor et al 
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ifTTl . therefore turning the 2D analysis into a 3D analysis. 

A further limitation of this technique is that we have assumed that the structure of the skin 
biopsy removed is representative of the entire tensile test sample. In reality, the structure, and 
therefore the properties of skin may vary considerably over a small area. The variation of the 
mean fibre orientation over the volume of the biopsy was quantified, however, the size of the 
biopsy is very small relative to the tensile test sample and we cannot infer that the variation 
would be negligible. At a minimum, future studies should excise a number of biopsies along 
the length of the test specimen to investigate the variation of the structure. 

While this technique has been described as 'automated', there are still a number of 'manual' 
steps that must first be performed. The first manual task is the histological staining, combined 
with the mounting of skin biopsies. The collagen detection process demands a high quality 
of histological staining for the method to be successful and therefore, care must be taken to 
follow standard procedures carefully. The second 'manual' task is the image acquisition phase. 
Modern 'slide scanners' automatically scan multiple slides at once meaning that tedious image 
acquisition techniques using a manual microscope are no longer necessary, however images 
must still be captured from the 'digital slide'. 

While this technique has provided quantitative structural data of human skin, it can, of 
course, only be applied in-vitro. Considering the effect that the mean orientation of collagen 
fibres has on the mechanical response of skin, the development of in-vivo methods for establish- 
ing the orientation of fibres is of the utmost importance. Advanced imaging techniques such as 
ultrasonic surface wave propagation may eventually provide real-time, in-vivo structural data. 

It should be noted that the nonlinear curve fitting technique rests only on measurements of 
Ai and CJn, and that A2, along with ki and k2 are obtained during the simultaneous optimisa- 
tion of Eq. ([To]) and Eq. ( pTj ). Ideally, a more complete analysis would include, compare, and 
contrast experimental data for A2 and/or A3. An extensive experimental data set would include 
planar biaxial tests with in-plane shear and separate through thickness shear tests ll9l lfT5lllfT8ll . 
however in the absence of these advanced testing protocols tensile tests coupled with a histolog- 
ical study of the collagen fibre alignment can be used for reasonable determination of material 
parameters [.14.1 . Of course non-uniqueness of optimal material parameters is an intrinsic prob- 
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lem in non-linear fitting. It is possible that the optimisation procedure finds a local minimum 
and assumes that this is the global minimum [25]. To ensure that the optimisation procedure 
is providing a unique set of material parameters a number of checks are available: The prop- 
erties of the Hessian matrix can be investigated at the optimum lfTT1lll2T1l . or alternatively, one 
can plot the objective function as a function of the varying material parameters. In this case 
the objective function was investigated and the plot (not reproduced) shows that our procedure 
calculates a global minimum and not merely a local minimum. 

In this study we have developed a simple automated process which can detect the orientation 
of collagen fibres. This technique can be easily implemented in MATLAB and can be adapted 
to detect other biological features, such as certain cells, leading to applications in diagnostics. 
We have applied this technique to skin biopsies and provided new quantitative data on the ori- 
entation of collagen fibres in the human dermis. So far, the availability of accurate structural 
data has lagged behind the progress of anisotropic constitutive modelling. Here we have pro- 
vided the structural data required to accurately make use of advances in constitutive modelling, 
and help fill the void of experimental data. The model parameters of the GOH model have been 
evaluated for skin using experimental data from the same skin samples. These sets of parame- 
ters will provide invaluable data for those wishing to model the anisotropic behaviour of skin. 
Finally, an FE simulation of a uniaxial tensile test on three separate human skin samples was 
performed which predicted the the response of these three samples well. We have illustrated 
that the Gasser-Ogden-Holzapfel model can successfully model the anisotropic behaviour of 
human skin and that it can be implemented in ABAQUS with ease. 
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6 Tables & Figures 



Table 1 : Local mean orientation of fibres and dispersion factor. Orientations are given with 
respect to the axis perpendicular to the axis of applied stretch. The three samples highlighted 
are those taken as illustrative examples for further analysis. (Note that the data given is axial 
data i.e. it represents undirected lines and does not distinguish between G and n + G [fT6l e.g. 
the orientation of ° and 180 ° are equivalent). 



Age 


Gender 


Location 


Orientation of 
Langer lines 

o 


Preferred 
Orientation 

o 


Dispersion factor K 


81 


Female 


5 
4 
6 
1 


0/180 
0/180 

45 
45 


174±3 
20±3 
38±7 
46±8 


0.1306±0.0054 
0.1439±0.0088 
0.1314±0.0054 
0.1675±0.0023 




^^^^1 


■_ 90 


88±8 1 


^0.1535±0.0059 




3 


135 


121±5 


0.1485±0.0026 


89 


Male 


5 


0/180 


178±4 


0.1462±0.0053 


^^^^ 


■ 0/180 ^ 


^1 0±5 


0.1456±0.0055 


6 

1 

2 


45 
45 
90 


61±4 

13±3 
89±5 


0.1289 ±0.0046 
0.1276±0.0054 
0.1009±0.0085 


3 


135 


118±7 


0.1602±0.0095 
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Table 2: Values obtained through curve-fitting for the parameters ^u, k\ and k2. Also displayed 
is R^, a measure of goodness of fit, and K and 7 obtained directly through histology. 



Sample 


MPa 


^1 


ki 


r 







K 


/?2 

% 


Parallel 


0.2014 


243.6 


0.1327 


41 


88 


0.1535 


99.54 


Perpendicular 


0.2014 


243.6 


0.1327 


41 





0.1456 


97.96 


Non- symmetric 


0.2014 


243.6 


0.1327 


41 


118 


0.1602 


94.40 
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(a) Original histology slide, (b) Binarised image after au- (c) Binarised image after 
Scale bar is 1mm. tomated thresholding. erosion step. 




(d) All identified collagen (e) Remaining bundles (f) Best fit ellipse about each 
bundles outlined in green. which meet area and fibre that meets the specified 

eccentricity criteria. criteria. 

Figure 1 : Images output from automated algorithm. 
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Figure 2: Histogram of collagen orientations. The two distinct peaks correspond to the pre- 
ferred orientation of the two fiber families. The angle, 7, is half the distance between the two 
peaks i.e. 7=41°. 
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Figure 3: Lattice structure of crossing collagen fibres with fibre dispersion taken into account. 
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Figure 4: Biopsies of skin samples for purpose of histological staining. Note that the biopsies 
have been sliced in three orthogonal planes. 
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Figure 5: Manual segmentation of collagen. Elongated fibres were marked by black lines, and 
their orientation was later measured manually. 
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(a) K-^ 0.0085 



(b) K = 0.25 




(c) K = 033 

Figure 6: Three dimensional representation of the orientation of collagen fibres. 
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Figure 7: Dimensions of test specimen (mm). 
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Figure 8: Location of tensile test samples shown in Table [T] (figure amended from llT9ll ). 
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Figure 9: Nonlinear curve fitting to obtain the constitutive parameters: ix and ^i are related to 
the early (infinitesimal) stress-strain part of the graph, see (a); k2, to the rest of the curve. 
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Figure 10: Comparison of sample parallel to Langer lines and perpendicular to Langer lines. 
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S, S11 

(Avg: 75%) 

^r +6.281 e+07 
+5.784e+07 
+5.287e+07 
+4.791 e+07 
+4.294e+07 
+3.797e+07 
+3.301e+07 
+2.804e+07 
+2.308e+07 
+1.811e+07 
+1.314e+07 
+8.176e+06 
+3.210e+06 



(a) 



(b) 



Figure 1 1 : Cauchy stress in Pa of sample strained by 50% (a) Parallel to the Langer lines (b) 
Perpendicular to the Langer lines. 
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S, S11 
(Avg: 75%) 

E+6.634e+07 
+6.082e+07 
+5.530e+07 
+4.978e+07 
+4.426e+07 
+3.875e+07 
+3.323e+07 
+2.771e+07 
+2.219e+07 
+1.667e+07 
+ 1.115e+07 
+5.632e+06 
+ 1.131e+05 




(a) 



S, S12 

(Avg: 75%) 
t- +2.035e+07 
I- +1.805e+07 

- +1.575e+07 

- +1.346e+07 

- +1.116e+07 

- +8.863e+06 

- +6.566e+06 

- +4.270e+06 

- +1.973e+06 

- -3.231e+05 

F-2.620e+06 
-4.916e+06 
-7.213e+06 




(b) 
Figure 12: Cauchy stress in Pa of a non- symmetric sample strained by 30% (a) CJn (b) Gu 
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Figure 13: Comparison between predicted model response and experimental force- 
displacement data for a non-symmetric sample. 
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